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TECHNICAL MEMORANDUM X -53868 


STUDIES IN SYSTEM SIMULATION 
SUMMARY 


The results of three mathematical studies of the feasibility of various 
methods of large scale digital system simulation are discussed in this report. 

A class of digital simulation formulae for large scale systems wherein 
numerical methods are applied directly to the system components is considered. 
Stability and error analyses are undertaken and function analytic techniques are 
employed to derive a class of stable, error preserving numerical methods for 
such simulation formulae. Techniques for implementing the component discre- 
tization process for systems containing both continuous and discrete components 
are indicated. 


The simulation of a dynamical system by a lower order system is con- 
sidered. The problem is formulated in terms of the Taylor series coefficients 
of the system impulse response from which lower order Fade approximations 
may be obtained by computationally appealing formulae . 

A data compression scheme for integer matrices wherein prime numbers 
are used for addressing is formulated. The technique allows a vector with 
integer valued entries to be characterized by a single positive number that is a 
linear function of the given vector and from which the column vector can be 
uniquely recovered by purely algebraic techniques. The procedure is computa- 
tionally limited only by the magnitude and storage accuracy required for the 
characterizing integer and is ideally suited for the characterization of topologi- 
cal matrices, which because of their inherent sparseness, are not affected by 
these limitations. 


INTRODUCTION 


The advent of ever larger systems over the past decade has necessitated 
significant changes in the methods and theory of system simulation. Systems 



with a handful of components and 10 or 20 degrees of freedom have given way to 
systems containing hundreds of components and having thousands of degrees of 
freedom. As such the semi -intuitive simulation procedures that once sufficed 
are giving way to formal algorithmic procedures designed to make optimal use 
of available computational capacity. Such procedures require that one make the 
maximum possible use of the connectivity structure inherent in a system, i.e. , 
isolating large systems into smaller subsystems, minimizing the complexity of 
such subsystems, etc. The studies discussed in this report are aimed to the 
further development of such techniques. 

This report is composed of three such studies of the theory underlying 
various techniques of large scale system simulation. Although formulated with 
computational feasibility in mind, these studies are mathematical, not computa- 
tional, in nature. The results have not been implemented nor have computer 
programs been written. Rather, mathematical studies of the feasibility of 
certain approaches to the large scale system analysis problem are undertaken. 
Existence theorems and series representations are used throughout the deriva- 
tion of the various results, though once formulated the various results can be 
implemented by purely algebraic techniques without iteration or root finding. 

All three studies are motivated by the realization that practical large 
scale simulation must deal with the components of a system as separate identities 
that are combined only after considerable analysis is undertaken at the compo- 
nent level. The first study, therefore, deals with the problems associated with 
applying numerical techniques directly to the system components, and a class 
of numerical methods that circumvent the stability problems inherent in such 
an approach are delineated. The second study is concerned with the possibility 
of approximating a subsystem by one of lower order prior to combining it with 
the remaining system components. Finally, the third study deals with a data 
compression technique for storing the large topological matrices encountered 
in large scale system simulation. 


DIGITAL SYSTEM SIMULATION 
BY MEANS OF COMPONENT DISCRETIZATION 


Introduction 

A differential system is typically composed of a collection of dynamical 
components characterized by (in general nonlinear) differential equations 
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(la) 


Dx. = f. (x.,I., t) 

1 L 1 1 


O. = h.(x. ,I.,t) 

i iii 


i = 1, 2, , . . , m 


(lb) 


together with a set of algebraic connection equations of the type 

K(I,0,U,W,t) = 0 (2) 

Here "D" is the differential operator, I = col(I.) and similarly for O and 
X while U and W are system input and output vectors respectively. Note 
that the equations induced by the common connection models (block diagram, 
linear graph, etc.) are always linear. It is, however, often convenient 1 to 
include the effects of nonlinear algebraic components in the connection ( rather 
than the component) equations; hence we allow for the possibility of nonlinear 
connection equations. 

Such a system is typically simulated by combining the component equa- 
tions with those for the connections to form a differential equation characterizing 
the entire system 2 ’ 3 

DX = F(X,U,t) (3a) 

W = H(X,U,t) , (3b) 

which may be simulated by standard numerical methods. Although the process 
of forming equations ( 3) from (1) and (2) is often complex, we may symboli- 
cally write 

F(X,U,t) ^ K F [f.](X,U,t) (4) 


and 


H ( X , U , t) =± K H [h.](X,U,t) 


(5) 


1. N. Prasad and J. Reiss, The Digital Simulation of Interconnected Systems 
( in preparation) , 

2 . Ibid . 

3. R. Saeks, State Equation Formulation for Multipart Networks, IEEE 
Transaction on Circuit Theory (to appear Feb. 1970). 
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where K and K are function valued (F) functions of a function variable 

r ri 

(f.) determined by the connection equations. The digital simulation of such a 


system thus takes the form 


A(E)X = F(X,U,t) = K F [f.](X,U,t) 

(6a) 

W = H(X, U, t) = K H [h. ](X,U,t) 

(6b) 


where A(E) is a function of the shift operator, E , which in some sense 
approximates the derivative, and F = K is a corresponding approximation of 

r 



An alternate approach to the simulation problem is to apply numerical 
methods directly to the component differential equations, ( 1) , yielding discrete 
approximations to the given differential components characterized by the differ- 
ence equations 

X. ( E)x. = f.(x.,I.,t) (7a) 

1 1 i i l 

i = 1, 2, , . . , m 

O. = h.(x.,I,t) (7b) 

iiii 

where A. ( E) and f. are approximations, as before, of D and f. . Now 

observing that the functions and , being determined entirely by the 

connection constraints, are unchanged by the component discretization, we 
obtain a discrete simulation of the overall system as 


X( E) X = K F ff.](X,U,t) 

(8a) 

W = K H (h.](X,U,t) 

(8b) 


where X(E) = diag [A.(E) ] . This approach to the simulation problem is 

termed digital simulation by component discretization and is the subject of 
the present work. 
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The component discretization concept is not new, having long been 
employed in the special case of linear time -invariant systems wherein one 
replaces the "Laplace transform" of a differential system by the "z-transform" 
of a discrete approximation. The difficulty with these approaches, and the 
basic problem to be considered in the sequel, is two-fold. First, a "good" 
approximation at the component level may not yield a "good" system approxi- 
mation, and secondly, stable component approximations may yield an unstable 
system approximation. This latter problem has been considered by Fowler [1] 
for the case of linear time -invariant systems composed of a single loop, but 
no general solution has been given. 

In the following sections the concepts of approximation error and stability 
of a numerical method are formalized, and necessary and sufficient conditions 
for a numerical method applied to the components of a system to yield a stable 
system equation are given. Finally, in the last section various techniques for 
implementing a digital simulation program by component discretization are 
indicated, In particular the simulation of systems containing both continuous 
and discrete components is considered and the possibility of using different 
numerical methods for different components is explored. 


Stability Analysis 

A numerical method is defined as the substitution of the difference 
equation 

\(E)x = f(x,I,t) (9) 

for the differential equation 

Dx = f(x, I, t) . (10) 

1. Definition. A numerical method is component stable if the difference 
equation, (9) , resulting from the substitutions 

A(E) — D 

and 

f — f 

is stable whenever the given differential equation, (10), is stable. 


5 



2. Definition. A family of numerical methods is system stable if each 
numerical method is component stable and moreover, the system difference 
equation 

X(E)X - K F [f.](X,U,t) 


resulting from the substitutions 
A.(E) — * D 

f. f. 

1 l 


i = 1 , 2, . . . ,m 


in the component differential equations is stable whenever the system differen- 
tial equation 

DX = K F [f.](X,U,t) 


is stable. 

The concept of component stability has been widely studied for both the cases 
of ordinary and partial differential equations. In the particular case of a linear 
time -invariant system with 

\(E) = E - 1 (11) 

necessary and sufficient conditions on the "approximation" f of f to assure 
component stability are known [ 2] . In the case of system stability, where one 
requires that numerical methods applied to the components of a stable system 
yield a stable difference equation characterizing the entire system, a " theory" 
is available only for single loop systems [1] . In the following, necessary and 
sufficient conditions for a numerical method to be system stable are derived, 
and the error inherent in the resulting approximation is analyzed. 

Initially we observe that if one allows unrestricted connections (except 
for the requirement that a differential equation, such as (6) , characterizing the 
entire system exist) , the function K may have arbitrarily large sentitivity. 

Thus any small variation in the f. will yield, for a sufficiently ill behaved 

K , to an unstable system difference equation. We therefore have: 
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3. Lemma. A necessary condition for a family of numerical methods 


ME) — > D 

1 

f. f. 

1 l 


i = 1, 2, . . . , m 


to be system stable is that 


f. = f. 

i i 


i = 1, 2, . , . , m 


Of course since system stability subsumes component stability we also have: 

4. Lemma. A necessary condition for a family of numerical methods 


\ (E) — ► D 
f. — f. 

i i 


i = 1, 2, . . . ,m 


to be system stable is that they are component stable . 

Consistent with the preceding lemmas it will be necessary to delineate 
those component stable numerical methods for which f. = f (that is a condition 

on the function X.(E) , which assures component stability when f. = f. is 

required) . Of course we must also assure that the operator A(E) approxi- 
mates the derivative operator if the solutions of the discretized system are to 
approximate those of the given differential system, and moreover, A^(E) must 

be chosen so as to assure the finite dimensionality of the discrete approximating 
system. For linear time -invariant systems the class of ?l(E) that satisfies 

these conditions is characterized by the following lemma: 

5. Lemma. Let an arbitrary linear time -invariant component be 
characterized by a stable differential equation. Then a necessary and sufficient 
condition on the numerical method 

\.(E) — D 
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to assure that the resulting discretized component will be characterized by a 
stable, finite dimensional difference equation whose solutions approximate those 
of the given component (differential equation) is that the function A. (w) , 

viewed as a complex valued function of a complex variable, satisfy the following 
conditions: 


a. A (w) is rational. 

b. A. ( w) approximates ln(w) . 

e. If |w | > 1 , A.(w) is analytic and ReA. ( w) > 0 . 

Proof: Condition a is the usual necessary and sufficient condition for a 
difference equation to be finite dimensional while condition b assures that A. (E) 
approximates D for if 

AHw) w ln(w) (12) 

then 

A.(E) » ln(E) = D (13) 

Here the validity of the operator approximation (13) follows from the complex 
function approximation by means of the spectral mapping theorem [3] and the 
operator equality of equation ( 13) is the usual semi-group representation of the 
shift operator through its infinitesimal generator, D [3] . Finally condition c 
is a necessary and sufficient condition for component stability. This may be 
demonstrated by representing the component differential equation by its "trans- 
fer function, " H^(p) , in which case the substitution 

A (E) - — * D (14) 

corresponds to the complex function substitution 

A. ( z) - p (15) 

where p is the "Laplace Transform" variable and z is the "z-Transform" 
variable. Clearly the difference equation resulting from the substitution ( 14) 
has "transfer function" 
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(16) 


H.(z) - H.[X.(z)] 


Now H.(p) is stable if, and only if, it is analytic in the region Rep > 0 while 
H. ( z) is stable if, and only if, it is analytic in the region |z| > 1 . If condition 
c is satisfied and the given differential equation is stable, EL ( z) is the composi- 
tion of a function taking the region I z I > 1 analytically to the region Rep > 0 
and a function analytic in this region. EL ( z) is therefore analytic in the region 

|z| > 1 and thus represents a stable difference equation, as required. Con- 
versely, if condition c is not satisfied, one can always find an EL(p) for which 

H.(z) does not represent a stable difference equation. The three conditions 

therefore combine to yield a component stable numerical method of the required 
type while the failure of any one implies that the numerical method is either 
unstable, infinite dimensional or does not approximate the given differential 
equation. 

It should be noted that the above lemma yields a condition that assures 
that every stable differential equation will be approximated by a stable difference 
equation. There are, of course, many values of function \.(E) which do not 

satisfy the conditions of the lemma yet still take some stable differential equa- 
tions to stable difference equations. For instance the forward difference formula 
takes stable differential equations with eigenvalues in a restricted region to 
stable difference equations [2] , but it does not take all stable differential 
equations to stable difference equations . Although the above proof of the lemma 
is valid only for the linear time -invariant case, we conjecture (but have not 
proven) that the result holds in general. 

The result of lemma 5, which gives a necessary and sufficient condition 
for a class of numerical methods to be component stable, also, upon combina- 
tion with lemmas 3 and 4, yields a necessary condition for system stability. 

In fact, this condition is necessary and sufficient. 

6. Theorem. Let the stable components of an arbitrary stable linear 
time -invariant system be discretized by the substitutions 


X(E) — D 
f. — - f. 

i l 


i = i, 2, . . . , m 


9 



Then a necessary and sufficient condition to assure that the resultant discretized 
system difference equation is stable, finite dimensional and has solutions that 
approximate those of the given system differential equation is that 

f. - f. i = 1, 2, . . . , m 

11 

and A.,(w) , viewed as a complex valued function of a complex variable, satis- 
fies the following three conditions. 

a. A. (w) is rational. 

l 

b. A.(w) approximates ln(w) . 

c. If |w| > 1 , A^(w) is analytic and ReA^(w) > 0 . 

Proof: The necessity of the conditions follows from the preceding lemmas. To 
demonstrate the sufficiency we must show that the discrete system difference 
equation resulting from the process of the theorem is system stable, finite 
dimensional and has solutions that approximate those of the given differential 
system. The latter two properties follow from the corresponding properties 
of the discretized components and the sufficiency of lemma 5. Finally for 
stability the sufficiency of lemma 5 assures that each component is stable while 
the discrete system difference equation induced by the component discretization 
is 


A( E) X - K_[f l(X,U,t) (17) 

-f 1 

where A(E) = diag[A.(E)] . Now equation ( 17) is just the discretization of the 

given system differential equation, viewed as a single component, by means of 
the numerical method 


A(E) — d 

(18) 

V'i' - w 

(19) 


Since each A. (E) satisfies the conditions of lemma 5 so does their direct 
i 

product, A(E); hence the numerical method of equations ( 18) and (19) is 
component stable (by lemma 5) and therefore takes the stable system differen- 
tial equation to a stable system difference equation. The specified collection 
of numerical methods is therefore system stable and the theorem is proven. 
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The theorem completely delineates the class of numerical methods that 
can be successfully employed for digital system simulation by component dis- 
cretization. Two questions, however, remain. First, what is the effect of the 
component discretization process on simulation error? Secondly, do any numer- 
ical methods satisfying the conditions of the theorem exist? In the former case 
we observe that the system difference equation obtained by the component dis- 
cretization process (17) is precisely the same difference equation that would 
have been obtained by applying the numerical method of equations (18) and (19) 
directly to the system differential equation. We therefore have: 

7. Corollary. The simulation error resulting from the process of 
theorem 6 is identical to that which would result upon applying the numerical 
method 


ME) — ~ D 

vv - w 

directly to the system differential equation. 

Finally it must be determined whether any numerical methods satisfying the 
conditions of the theorem exist. This is indeed the case and, in fact, both 
trapezoidal 

\.(E) = (E + 1) _1 (E - 1) (20) 

and backwards 

X.(E) = E" 1 (E - 1) (21) 

integrations satisfy the required conditions as does the second order scheme 

\ (E) = ( 2E 2 ) _1 ( 3E 2 + 1) (22) 

Of course numerical methods of arbitrarily high order that satisfy the required 
conditions can be obtained by standard approximation techniques. It is note- 
worthy that the conditions of the theorem are never satisfied when X. is a 

polynomial in E , as is the case for Simpson's rule and most standard integra- 
tion techniques. 
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Digital Simulation 


A digital simulation program employing component discretization has 
two distinguishing characteristics. First, since one discretizes continuous 
components as a first step in such a program, components that are discretely 
specified may be included in the system simply by skipping the initial discreti- 
zation step. Secondly, since each component is discretized independently, 
different approximations [A.(E)]. may be used for different components. One 

may therefore use "better" numerical methods for simulation of those components 
that have small time constants and/or whose behavior is highly sensitive than 
for the remaining components. 

A component discretization simulation program might allow for four 
classes of dynamical components along with the usual memoryless and connection 
components. These include continuous finite dimensional components specified 
by an ordinary differential equation, continuous ( possibly infinite dimensional) 
components specified by a convolution integral, discrete components specified 
by either a difference equation or a discrete convolution, and finally components 
specified by a computer program, such as might arise if data analysis techniques 
are used to estimate the dynamics of an unknown device. Clearly if such a sys- 
tem is to be digitally simulated, one must use component discretization since 
some of the components are specified discretely and others (characterized by 
convolution integrals) cannot be stored on a computer in continuous form even 
though they are specified continuously. 

In implementing such a program, one would input the discrete components 
as given and the continuous components together with a specified numerical 
method, i.e. X.(E) , satisfying the conditions of theorem 6. Now as a first 

step the program discretizes the continuous components, The finite dimensional 
components characterized by differential equation are discretized by making 
the substitution 

X.(E) — * D (23) 

and then converting the resultant higher order difference equation to a first 
order difference equation by standard methods. On the other hand those com- 
ponents specified by convolution integrals are discretized by converting them 
to a discrete convolution with weighting coefficients determined by the numeri- 
cal method. Once all of the components have been discretized, the remainder 
of the simulation may be carried out by standard discrete system methods. 
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It is interesting to note that the process of discretizing a continuous 
component by the numerical methods of theorem 6 is essentially a matter of 
rearranging and indexing the matrices characterizing the given continuous 
component and involves no "real computation. " Additional computation, corres- 
ponding to the degree of the numerical method, is, however, reflected in the 
readout process wherein one is dealing with a difference equation of higher 
order than the given differential equation. For this reason no computational 
savings is gained by carrying out the discretization at the component level 
since the readout process is still carried out at the system level. This is not 

the case for those numerical methods with f. =£ f. , in which case the calcula- 

11 

tion of f. is usually more easily carried out at the (essentially decoupled) 

component level than at the system level, but such processes do not have the 
system stability characteristics required to successfully carry out the compo- 
nent discretization process. 


Conclusions 


Although we have carried out a rather long and tedious derivation, the 
results of the theory are readily applicable. As long as one employs integra- 
tion techniques of the type indicated by theorem 6, system stability is assured 
independently of the connections or type of components in a system. In fact, 
even if one integrates with too large a step size, the resultant solution, though 
inaccurate, will tend toward zero rather than becoming unstable. 

It is interesting to note that in a number of computer aided network 
analysis programs wherein component discretization is used, it has been found 
experimentally that trapezoidal integration is stable while Simpson's rule and 
other polynomial integration routines may be unstable. This is, of course, 
verified by our theory, which also yields a means for obtaining higher order 
stable integration routines. Another area wherein component discretization is 
employed is in the construction of a Digital Difference Analyzer . In such a 
device one immediately replaces integrators in the program with a discrete 
approximation independently of the remainder of the program and is therefore 
using component discretization, if implicitly. As in the case of network analysis, 
it has been found experimentally that trapezoidal integration is stable, as indi- 
cated by the present theory. 
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PADE APPROXIMATION AND 
THE DEGREE REDUCTION PROBLEM 


Introduction 

In system analysis it is commonly desired to replace a given system by 
one that has a lower degree ( order) but whose external behavior is similar or 
identical to the given system. Such problems are termed "degree reduction 
problems" and are the subject of the present work. 

Possibly the most common manifestation of the degree reduction problem 
in system analysis is the problem of eliminating uncontrollable and/or unobserv- 
able modes from a system to be simulated. That is the removal of internal 
responses that have no effect on the overall external system behavior. Although 
such modes represent unlikely situations for a single component, they occur 
commonly in interconnected systems wherein the effect of a mode in one com- 
ponent is canceled by an equal and opposite effect in another component. For 
instance in an electric network two parallel capacitors exhibit one rather than 
two independent modes. 

A second class of degree reduction problems encountered in system 
analysis is concerned with modes that have a small but nonzero effect on the 
overall system behavior. Such characteristics are often encountered in control 
systems wherein a mode is nominally unobservable but, because of component 
variations and/or simulation error, appears in the output with a small residue. 
For instance a plant mode that has nominally been canceled by compensation 
techniques may in fact appear with a small residue because of variations of 
the components from their nominal values. In such a case the computational 
saving achieved by neglecting such a mode may justify the error induced into 
the simulation. 

Finally practical considerations, such as the complexity of the system 
or the amount of memory required for simulation, may force one to approximate 
a given system by one of lower degree even at the cost of considerable simula- 
tion error . We thus have a third class of degree reduction problems — approxi- 
mation. 


The study of these three classes of degree reduction problems in a 
context amenable to digital system simulation is the purpose of this work. 
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Since the theory is essentially a reformulation of the realizability theory of 
Youla [ 4, 5] , Ho [6, 7] , and Kalman 4 , we state the main results without proof, 
referring to the recent text of Kalman, Falb, and Arbib [8] . In the following 
a "universal system specification" is formulated wherein a system is charac- 
terized by an infinite sequence of matrices . This sequence, which has 

only a finite number of independent terms and is thus computationally feasible, 
can be obtained by inspection from the system impulse response, transfer 
function, and state equations or directly from measured data and thus serves 
as a natural intermediary between the various system specifications. Once 
such a sequence has been constructed (starting with any of the usual system 
characterizations or measured data) , minimal state equations for the system 
that solve the various degree reductions problems are obtained. Since the 
entire procedure is algebraic, the resulting algorithms are ideally suited for 
digital system simulation, no approximation or iteration steps being required. 


System Specification 

The most common methods for characterizing a linear time -invariant 
(finite dimensional) dynamical system are the state equation 

DX = FX + GU (24a) 

Y = HX (24b) 

the system impulse response, K(v) , such that 

OO 

Y(t) = f K(t - q)U(q) dq (25) 


and the transfer function, T(p) , such that 

Y(p) = T ( p) U ( p) (26) 

where Y(p) and U(p) are the "Laplace Transforms" of Y(t) and U(t) 
respectively. Now it is well known that these characterizations are equivalent 
yet the interrelationship between the various characterizations are quite com- 
plex. If one, however, expands K(v) and T(p) in appropriate series, a 
commonality between the three characterizations can be found. Indeed this 
commonality corresponds to readily measurable parameters of the system. 


4. B, L. Ho and R. E. Kalman, the Realization of Constant Input -Output Maps, 
SIAM Journal on Control (to appear) . 
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Consider a system impulse response matrix, K(v) , which if it 
represents a finite dimensional dynamical system, has a Taylor series expansion 

o° 

K(v) = (27) 

1=1 1 

about the point v = 0. Here the ith coefficient matrix, A. , is the (i-l) st 

derivative of K(v) evaluated at v = 0 and the coefficient sequence completely 
characterizes the system. Although as defined, the A, sequence appears to 

be infinite, only a finite number of the terms are independent; hence the values 

of A. form a computationally feasible characterization for a finite dimensional 

dynamical system. In fact, we have the following fundamental result, which 
completely characterizes the finite dimensionality of the A. sequence. 

1. Theorem. Let K(v) be the impulse response of a finite dimensional 
system. Then there exists a set of real constants, b 1 ,b 2 , . . . .,.b , such that 


A. = 

l 


- X b, A. . 4 

, Lj . k k+i-n-l 
k=l 


i = n+1, n+2, . . . 


Moreover, the smallest n for which this is true is the dimension of the mini- 
mal system ( i.e. number of state variables) that realizes K(v) exactly, and 
the minimal polynomial of such a system is 


\{z) 


n , n-1 , n-2 , 

= z + bjz + b 2 z + . . . + b 


A short proof of this fundamental theorem is given in Reference 8. 
Consistent with the theorem, a knowledge of the first n A.'s together with the 
n b^'s is sufficient to completely characterize the system. Equivalently the 
values of the first 2n A.’s completely characterize the system since the values 
of b.'s can be calculated from these by means of the equality of the theorem. 
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Although the values of A. have been derived in terms of the system 

impulse response, they have an equally natural interpretation in terms of the 
system transfer function. In this case one takes a Laurant expansion of T(p) 
about infinity and obtains the series expansion 

co 

T(p) = 2 a /p l (28) 

i=l 

The fact that the Laurant coefficients are indeed the same as the Taylor coef- 
ficient of the impulse, response follows immediately upon an application of the 
initial value theorem, Of course theorem 1 holds whether the A^ sequence 

is obtained from the impulse response or the transfer function. 

Clearly the A^ sequence serves as a natural intermediary between the 

frequency and time domains, and in fact completely characterizes the properties 
of both. In fact the A. sequence is also naturally related to the state charac- 
terization of the system. To see this we observe that for any state equation that 
realizes a given (external) system 

K(v) = He FV G , (29) 

which upon differentiating and evaluating at zero yields 

A. = HF 1_1 G (30) 

The solution of the converse problem of identifying a state equation, such 
as (24) , or equivalently the three matrices H, F, and G, given a sequence of 

A. , is by no means obvious. Fortunately it has a computationally appealing 

solution, though one that requires a rather tedious derivation. This result, 
due to Youla [4] and Ho [6, 7] is predicated on the properties of a Hankel 
matrix associated with the A. sequence. We therefore let J be the block 
symmetric matrix. 
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Ai A 2 A 3 . . , 

A 

n 

A 2 A3A4 . . • 

A < 
n+i 

A3A4A5 

• 

• 

• 

• 

• • 

• 

• 

A A 

A 

n n+i 

2n- 


where n is any integer such that the condition of theorem i holds. Similarly 
we define J' to be the same matrix with the subscripts shifted up by one, 
i.e. , the 1-1 entry in J' is A 2 and the n-n entry is A . With this mode 

of specifying the A sequence, a minimal state equation realizing the system 

characterized by the A. sequence exactly may be obtained by purely algebraic 
manipulation. 

2 . Theorem. Let a system have Hankel matrix, J , and let P and M 
be arbitrary nonsingular matrices such that 


PJM 



0 I 0 


Then 


DX = FX + GU 
Y = HX 

is a state equation of minimal dimension realizing the given system exactly if 
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H = 

[i ! c 

f| JM 

fl 1 

n 


L P» 

J 

_ 0 _ 


Here I is the n by n unit matrix, n is the rank of J , I and I are 
n J r P 

r and p dimensional unit matrices, r is the number of system inputs, p is 

the number of outputs, and in the partitioned matrices, 0 is the zero matrix of 

conformable dimension. 

Although the proof of the theorem is quite involved (see for instance 
Reference 8) its application is quite straightforward. Operationally J is formed 
from the sequence and is diagonalized by P and M . These matrices are 

then used to calculate F, G, and H . Although it is necessary to diagonalize 
J, P and M need not be in any way related; hence the diagonalization may be 
carried out by independent row and column reduction processes, which yield 
the diagonalization in a fixed number of steps without iteration or root determi- 
nation. 


If one knows a priori a minimal set of b^ , as in the next section, such 

that theorem 1 is satisfied, an even simpler state realization that eliminates 
the diagonalization is possible. 

3. Theorem. Let a system have Hankel matrix, J , and let b^ , 

i=l, 2 n be a minimal set of real numbers such that the equality of theorem 

1 holds. Then 

DX = FX + GU 

Y = HX 

is a state equation of minimal dimension realizing the given system exactly if 


19 




Consistent with the preceding development the A sequence, or more 

accurately its first 2n terms, serves naturally as a universal system specifi- 
cation for linear time -invariant (finite dimensional) dynamical systems. The 
sequence can be obtained from either the time, frequency or state representa- 
tion or alternatively from direct measurements of the system. Conversely any 
of the three system representations can be obtained from the A. sequence by 

purely algebraic manipulation. Moreover, the close relationship between the 
A. sequence and measurable system parameters (the impulse response, 

frequency response, etc. ) renders it an ideal medium in which to carry out 
degree reduction and/or system approximation. 

Degree Reduction Problems 

With the formulation of the A. sequence and its relationship to the 
common system models we may proceed to attack the various degree reduction 
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problems. In fact, much of the preceding theory was developed with the solution 
of the first degree reduction problem in mind (i.e. , the problem of eliminating 
uncontrollable and/or unobservable states from a state equation) . This is 
achieved by starting with an arbitrary state equation representing a given sys- 
tem, such as 

DX = FX + GU (32a) 

Y = HX , (32b) 

and constructing (the first 2n members of) its sequence as per equation 

(30) . Now J and J' are constructed from this sequence and used in theorem 
2 to obtain a minimal state equation having the same input-output behavior as 
the given equations. 

DZ = FZ + GU (33a) 

Y = HZ . (33b) 

As such, one can construct a pair of state equations that have the same external 
behavior as the given equations but without the unobservable and uncontrollable 
states; hence the first degree reduction problem is solved. In fact, the solution 
is entirely algebraic. 

Unlike the first degree reduction problem, the second and third problems 
require that one find a state equation whose impulse response approximates that 
of a given higher order system. One must therefore choose a mode of approxi- 
mation before a solution can be obtained. There are, of course, many possible 
approaches to the approximation problem; for instance, L* , L 2 , Chebychev, 
etc. , each of which has merits in various contexts. In the present problem a 
Pade approach [9] will be taken, primarily because of the computational simplic- 
ity inherent in such an approach; that is, an impulse response K(v) is said to 
be a kth order Pade approximation of an impulse response K(v) if the first 
k Taylor series coefficients of K(v) coincide with those of K(v) . In terms 
of the A. sequence we therefore require that 

A. = A. i = 1,2, . . . ,k . (34) 

Consistent with the close relationship between the A. sequence and the various 
system models, such a mode of approximation is computationally appealing 
(if not of great physical relevance) . 
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From theorem 1 it follows that given any arbitrary set of matrices 


A. i = 1, 2, . . . , n 

i 

and any set of n real numbers 

b. i = 1, 2, . . . , n 

l 


sequence defined by equation ( 35) for 


A. 

i 


n __ 

E b i A , • < 

. u . k k+i-n-1 
k=l 


(35) 


(36) 


(37) 


for i >n has an n dimensional state equation realization (for instance as 
obtained by theorem 3) . Moreover, the first n Taylor series coefficients of 

the corresponding impulse response are the specified values of A^ , and the 

minimal polynomial of the resulting realization is 


, , , n n-i n-2 , 

A(z) = z + b 1 z +b 2 z + . . . +b 


Clearly if one is given an arbitrary A. sequence, we can construct an nth 

order system, which is an nth order Pad& approximation to the given A. 
sequence, simply by letting 


A. = A. 
i i 


i = 1, 2, , . . , n 


and choosing the remainder of the A^ sequence by use of equation ( 37) for 

some arbitrarily specified set of b. . In fact, since the values of b. are 

arbitrary, one can completely specify the system dynamics and still achieve 
the required approximation. Formally we have: 

4. Theorem. Let an arbitrary system be characterized by an A. 
sequence . Then for any set of coefficients 

b. i = 1, 2, . . . , n 

i 
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the sequence 



i = 1,2,. 


,n 


A. 

l 


n 

■ T. b A 

, u . k k+i-n-1 
k=l 


i = n+1 , n+2, . . 


characterizes an n-dimensional system having minimal polynomial 

^ , . n n-1 , 

M z) = z + b,z + . . . + b 
1 n 

which approximates the given system to the nth degree in the Pade sense. 

Consistent with theorem 4 the second and third degree reduction problems are 
solvable by algebraic means if one employs the Padd mode of approximation. 

It should be noted that the Padd mode of approximations assures that the approxi- 
mation of the impulse response will be "good" (near v = 0) while stability 
assures that both the actual and approximate impulse responses will tend to 
zero for large v ; hence the Padd mode of approximation is quite reasonable . 

Of course the quality of the approximation can be improved if one makes an 
appropriate choice of the characteristic polynomial A( z) , which is arbitrary 
except for the requirement that it be stable. One approach that further improves 

the approximation near v = 0 is to choose the b^'s so as to increase the order 

of the Pade" approximation. That is the values of b are chosen so that 

A. = A. n < i < m , ( 40) 

ii 

thereby obtaining an mth order approximation with an nth order system. The 
circumstances under which this can be done have been studied 5 and will not be 
delineated here. If appropriate rank conditions on J are satisfied, it is, 
however, possible to achieve a Padd' approximation with order as high as 2n . 

An alternative approach, which allows one to control the approximation 
when v is neither large nor small, is to choose X( z) so as to approximate the 
characteristic function of the given system. In this way one can guarantee that 
the dynamics of the approximate system (i.e. , oscillitory frequencies, decay 


5. B. L. Ho and R. E, Kalman, op. cit. 
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rates, etc. ) are similar to those of the given system. In particular for the 
second degree reduction problem where one simply desires to eliminate negli- 
gible but nonzero modes if A(z) is taken as the factor of the characteristic 
function of the given system corresponding to the nontrivial modes, then the 
dynamics of the approximate system will be the same as those- of the given 
system except for the deletion of the negligible modes (i.e. those with small 
residues) and slight variations of the residues to the remaining modes to com- 
pensate for the effects of the deleted poles at v = 0 . The use of Pad£ approxi- 
mations of a system impulse response thus yields a complete solution for the 
second degree reduction problem and a solution to the third degree reduction 
problem in those circumstances wherein the Pade’ criterion is relevant. More- 
over, in both cases the procedure is completely algebraic and computationally 
appealing. 


Conclusions 


Although the derivation of the preceding results is quite complex, the 
results lend readily to the development of computational algorithms for the 
solution of the degree reduction problem. The required steps to implement 
such an algorithm are as follows. 


1. First Degree Reduction Problem. Given an arbitrary state equation 
representing a given system, characterized by the three matrices F, G, and 

H we first form the sequence for the system via A^ = HF i *G . Now these 

values of A. are used in theorem 2 to form new minimal state equations, 

characterized by the matrices F, G, and H. Since this equation is minimal, 
all uncontrollable and unobservable modes have been removed while, according 
to the theory, the system (as observed externally) is unchanged. 

2. Second Degree Reduction Problem . Given an arbitrary state equation 
characterized by matrices F, G, and H , which has characteristic polynomial 
A( z) , let us assume that A( z) = Ai(z) A 2 ( z) where Aj(z) corresponds to n 
significant modes and A 2 (z) corresponds to m negligible modes. Now the 
second degree reduction problem is to find an nth order system that is the 
same as the given system except for the elimination of the negligible modes. 

First one forms the A. sequence via A. = HF* -G and lets A. = A. for i 

i i l i 

less than or equal to n . Now the values of b^ are taken as the coefficients 
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of X t (z) and the remainder of the A sequence is formed via equation (37) . 

Finally either theorem 2 or_3 is_used to_find the approximate nth order system 
characterized by matrices F, G, and H . 

3. Third Degree Reduction Problem. Given an arbitrary kth order 
system, characterized by matrices F, G, and H and characteristic polynomial 
A(z) , the third degree reduction problem is to find an nth order system that 
approximates the given system in some sense. Initially we choose T( z) to be 
an nth order polynomial that in some sense approximates A( z) . Now the 

sequence, A. = HF* *G , is formed and used to define A. = A. for i less 

l __ i i 

than or equal to n . Finally the remainder of the A^ sequence is obtained by 
equation (37) , using the coefficients of \(z) for the b^'s , and the state equa- 
tions for the approximating system, with matrices F, G, and H, are obtained 
either by theorem 2 or 3. The key to the third degree reduction problem is the 
choice of the approximation used to choose A( z) . A decision on this can only 
be made in the context of the application. One possible choice, however, is to 
use a Pade approximation about zero. Since this would assure that the low 
frequency characteristics of the given system were preserved by the approxi- 
mate while the Pad£ approximation of the A^ sequence assures that the high 

frequency characteristics are preserved, one would expect a reasonably close 
approximation. 

LINEAR DATA COMPRESSION FOR TOPOLOGICAL MATRICES 


Introduction 

The topological, or connection, information in a large scale system is 
often conveniently characterized by a connection matrix of some type. Such 
matrices typically have integer entries corresponding (in some sense) to the 
existence of a connection between two components and zeros in the remaining 
entries. Since in a large scale system a given component is typically connected 
only to two or three others, these matrices are generally sparse, though this 
sparseness may have no discernible pattern. With the ever increasing size of 
modern systems, such matrices may have tens of thousands of entries; hence 
some form of data compression is necessary if they are to be effectively 
employed in the digital simulation of such systems. 
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There are two basic approaches to the topological data compression 
problem. One is based on the fact that memory words are used inefficiently by 
topological matrices ( since many of the entries are zero and need not be stored) 
and the other based on the fact that available word length is not used efficiently, 
since the entries are usually taken from an alphabet of small integers (including 
zero) which do not require a full word length for storage. A number of data 
compression techniques based on the first approach are commonly used. Typi- 
cally these store only the nonzero entries together with an address and are often 
predicated on an a priori knowledge that the nonzero entries are arranged in 
some specified pattern (such as one nonzero entry per column, etc.) . It is 
a data compression scheme of the second type that is considered here wherein 
a single element from a large alphabet characterizes a vector of elements taken 
from a small alphabet. As such, one stores a large integer rather than a collec- 
tion of small ones, hence making more efficient use of available word length. 

In the theory neither sparseness nor "smallness" of the matrix entries is 
required, but computational considerations demand both; hence the technique is 
well suited for topological matrices but not generally applicable. 

In the following sections the data compression scheme is developed 
without regard to computational considerations (i.e. sparseness, magnitude of 
the entries, etc. ) , first for vectors with nonnegative integer entries, then for 
vectors with arbitrary integer entries and finally for matrices with arbitrary 
integer entries. In all cases the scheme is shown to be linear (either in the 
sense of a group or a semi -group) ; hence linear operations on the matrices or 
vectors may be implemented by carrying out corresponding operations on their 
compressed characterizations. Finally algebraic recovery algorithms are 
formulated and an error analysis is carried out wherein the limitations that 
must be imposed on the technique to assure accurate recovery of the given 
matrix via computational methods are delineated. 


Representation Theory 

Initially we will consider a column vector of nonnegative integers. 

The length of such a vector need not be restricted and hence these vectors may 
be taken as infinite if one requires, in addition, that they have only a finite 
number of nonzero entries (i.e. they are zero almost everywhere) . Such a 
vector is then a member of the space 


V 


QO 


2 N ° 

i=l 


(41) 
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(N° the nonnegative integers) which forms a semi -group under componentwise 
addition. A vector a in V is denoted by a - (a 1? a 2 , a 3 , . „ . ) .= (a^) . Also 

let p. be the ith prime number. That is p 4 = 2, p 2 ■= 3, p 3 = 5, p 4 = 7, p 5 = 11, 

etc. (An ordered list of prime numbers is given in Reference 10.) 

Now given any vector a in V , let 0(a) be the positive integer 

a. 

0(a) = 5 p. 1 (42) 

1=1 1 

Note that since a in V is zero almost everywhere all but a finite number of 
a. 

the factors p^ are unity; hence for all practical purposes the product of 

equation (42) is finite and can be readily computed. The possibility of identi- 
fying a single integer with a column vector of integers is certainly not surprising 
and can be done by any number of formulae. What is possibly more surprising 
in this case, however, is that the function 0(a) is a linear function on the 
semi-group V and moreover an isomorphism. As such, a can be recovered 
uniquely from 0(a) . To this end we have: 

1. Theorem. The function 

0 : V — N 1 

a — 0(a) 

defined by equation (42) is a semi-group isomorphism. (Here N 1 is the multi- 
plicative semi-group of positive integers) . 

Proof: To prove the theorem we must show that the function, 0 , is one to 
one, onto and linear. In the first case if 0(a) = 0(b) , then 

a. oo b. 

n p. 1 = 0(a) = 0(b) = n p. 1 (43) 

i=l 1 i=l 1 
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a. b. 

Now since only a finite number of the factors p. and p are not unity, the 
unique prime factorization theorem [10] applies to equation (43) ; hence a^ = b^ 
for all values of i showing that a = b , as required. Therefore 0 is one to one. 

To verify that 0 is onto N 1 , let n be an arbitrary positive integer 
with prime factorization. 


n 


m a i 
n P. 
i=l 1 


a. 

°° 1 

n p. 

i=i 1 


(44) 


where on the right side of equation ( 44) , a. is taken as zero if p. is not a 
prime factor of n . Clearly 


n = 0(a) = 0( (a.) ) 


and 0 is therefore onto. 


Finally if a and b are in V , 


(a,+b.) 
00 11 
0(a + b) = n p. 

i=l 1 



hence 0 is linear as required. 


(45) 


0(a)0(b) 


(46) 


Note that repeated application of the formula 0(a+a) = 0(a) 0(a) yields the 
scalar multiplication formula 

0(ka) = 0 ( a) k (47) 

for any nonnegative integer k . The theorem assures, at least in theory, that 
a vector a in V can be recovered uniquely from 0(a) , but no computational 
algorithm for the recovery process is indicated by the proof. In fact several 
computationally appealing algorithms are possible and will be formulated in 
the following section. 

The preceding theory can be extended to the case of vectors with entries 
composed of arbitrary integers if one adopts a rational, rather than integer, 
representation. In this case one is dealing with vectors from the space 
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u 


N 


(48) 


_ V 
- U 
i=l 

of infinite vectors having a finite number of nonzero integer entries. Clearly 
U is a group under componentwise vector addition. 

For an arbitrary vector a in U let 

a = a + - a - (49) 

be a decomposition of a into the difference of two vectors from V and define 
0(a) as 


0(a) = 0(a + )/0(a“) (50) 

where 0(a + ) and 0(a~) are defined as per equation (43) for a + and a - in 
V . This is well defined for if a_ + and a~ represent the decomposition of a 
wherein the positive entries are taken for a + and the negative entries for a“ , 
any other decomposition has the form a + = a + + k and a“ = a - + k where k 
is in V ; hence 


0(a + )/0(a“) = 0(a + + k) /0 (a - + k) = 0(a_ + )0(k)/0(a_ )0(k) 

= 0(a+)/0(a“) , (51) 

showing that 0(a) is independent of the decomposition employed. Also 0(a) 
defined as per equation (50) is an extension of 0 (a) as defined in equation 
(42). That is if a is in V, then the two definitions of 0(a) coincide (by 
taking a - = 0) . An argument similar to that of theorem 1 will yield the 
following result wherein it is shown that 0 is a group isomorphism when 
defined on U ; hence the unique recovery and linearity properties obtained 
for 0 when operating on V also hold for its extension. 

2. Theorem. The function 

0 : U — Q + 

a — * 0(a) 

defined by equation (50) is a group isomorphism. 
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Although we have formulated the preceding theorem in terms of rational 
numbers, in practice one would probably store the numerator and denominator 
integers separately, in which case it would only be necessary to recover a 
vector in V from an integer and one would not have to deal with rational 
numbers computationally. To justify such an approach, however, it is necessary 
to show that a can be recovered uniquely from 0(a) independently of the 
integers used to represent 0 (a) . Assuming an algorithm is available with 
which to calculate 0 _1 (n) for any integer ( see the following section) , it suf- 
fices here to show that 

a = 0 -1 ( n) - 0~ 1 (m) (52) 

is independent of the integers n and m used to represent 0(a) = n/m . To do 
this let n and m be the unique pair [10] of relatively prime integers such that 
0(a) = n/m . Then any other such representation, such as 0(a) - p/q , satis- 
fies the equalities 

p = kn ( 53) 

and 

q .= km ( 54) 

with k in N 1 ; hence 

0 -1 ( p) - 0 -1 (q) = 0 _1 (kn) - 0 -1 (km) 

= 0 -1 (k) + 0 -1 (n) - 0 -1 (k) - 0 _1 (m) (55) 

= 0 -1 (n) - 0 -1 (m) 

The vector a recovered is therefore independent of the particular integers 
p/q used to represent 0(a) ; hence it is reasonable to store two separate 
integers from which common factors need not be removed before recovering a . 

The preceding vector representations can be extended to form a matrix 
representation in a number of ways. Possibly the most convenient is to simply 
"unfold" a p by r matrix to form a pr vector and apply the preceding repre- 
sentations directly. This, unfortunately, is not algebraically well behaved and 
results in extremely large integer representations, which, as will be shown in 
the sequel , is the primary limitation on the application of the scheme . We 
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therefore prefer to apply the representation separately to each column vector 
of a matrix separately, thereby yielding a representation for a matrix as a row 
vector of column representations. That is for a matrix, M , we have 

0(M) - row[0(m 1 ) , 0(m 2 ) 0( m r ) ] = row0(m ) (56) 

where m. is the ith column vector of the matrix M . So defined, 9 is a 
1 

natural extension of the preceding vector representation since if M is composed 
of a single column vector, the single representing rational number of equation 
( 56) is the same as that of equation ( 50) . Also the linearity of 0 is preserved 
via the formula 

9 ( M+N) = row.[0(m )'0(n )] (57) 


and the one to one, onto properties of 9 follow from the corresponding proper- 
ties for the individual column vectors. As before there is no limitation on the 
size of the matrices employed; hence we deal with infinite matrices that are 
zero almost everywhere, this group being denoted by W . Unlike the previous 
development, our representation is a vector of positive rational numbers, rather 
than a single number. Since our matrices are zero almost everywhere, the 
representing vectors are one almost everywhere. We therefore denote by S+ 
the space of infinite vectors with positive rational entries, which are one almost 
everywhere. Clearly this space is a group under componentwise multiplication 
for which we have the following representation theorem. 

3. Theorem. The function 

0: W — S + 

M — 0(M) 

defined by equation (56) is a group isomorphism. 

Since 9 is linear, one can add matrices by carrying out a component- 
wise multiplication on the entries of their representations. Unfortunately a 
similar formula for matrix multiplication does not exist. We can, however, 
multiply the representation of a matrix by an unrepresented matrix to obtain 
the representation of its product. A little algebra will reveal that for matrices 
M and N 
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0 ( MN) = row 


(58) 


" r 

n e(m.)n,j 

J=i j J1 


hence matrix multiplication can be done in a sense . 

Finally we note that one can carry out column operations directly on 
0 ( M) for column interchange; multiplication of a column by a scalar and the 
addition of two columns can all be carried out directly in terms of the 0(m.) . 


Recovery Algorithms 

In the preceding development we have shown that, in theory, a can be 
uniquely recovered from 0(a) . It, however, remains to obtain a computa- 
tionally feasible algorithm for calculating a given 0(a) , Clearly it suffices to 
consider the case where 0(a) is in N 1 and a is in V since the representa- 
tions for a in U and a in W are just collections of representation of the 
former type. Our fundamental lemma in this respect is the following wherein 
the square bracket notation is used to denote the least integer operation [10] . 

4. Lemma. For any integer n and prime number p 

[l - n/p k + [n/p k ]] 

is one if p divides n and zero otherwise, 
ic ic 

Proof: If p divides n , the n/p is an integer; hence 

n/p k = [n/p k ] (59) 

and 


[l - n/p k + [ n/p k ] ] = [1] = 1 (60) 

On the other hand if p does not divide n , 

n/p k - [n/p k ] (61) 
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is a fraction between zero and one (exclusive) ; hence one minus this fraction is 

also such a fraction and therefore its greatest integer is zero as required. In 

[ 1c k 

1 - n/p + [ n/p ] J is one if n/ p 

is an integer and zero if it is a fraction and is therefore readily evaluated com- 
putationally. Recognizing that the number 0(a) (for a in V) has a factor 
a. 

p. 1 , the lemma leads immediately to the following recovery theorem. 

5. Theorem. 


1 - 0 ( a) / p. + 


0(a)/p. k 


if and only if k s a- . 

Operationally there are a number of ways in which theorem 5 can be 
applied. If one is dealing with possibly large integers, a. , the most efficient 

t 1 k T kll 

i - 0(a) /p^ + 0(a)/p - i 

(i.e. the largest k such that .0(a)/p‘ is an integer) . Assuming that the 
integers a. are known to be bounded by 2 m , this search can be done in m 
steps, as follows. 


,m 


6. Corollary. Let a in V be such that a. < 2 , Then a^ can be 

calculated from 0(a) in no more than m steps by means of the following 
algorithm . 


a. 

Let d = 2 m 

1 

3 

let j 

= 1 . 




b. 

If 0(a)/p. d 

is 

an integer 

a. 

i 

2= d, 

go to c . 


If 0(a)/p. d 

is 

not an 

integer 

a. < 

i 

d , go to d. 

c. 

Let d = d + 

m-j-l 
2 J ; 

let 

j = 

j + * 

; go to b. 

d. 

Let d = d - 

m-j-l 
2 J ; 

let 

j = 

j + 1 

; go to b. 
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Clearly after m steps of the above algorithm the value of a. has been "trapped" 
and the recovery is complete. 

A second approach to the recovery process is well suited for the case 
when a^ is known to be small, as for instance in the storage of topological 

matrices. In this case a^ can be calculated by the following formula, which 

converges in a. + 1 steps. 

7. Corollary. 

1 - 0(a)/p k + j^0(a)/p. k 



Here if a term is zero, all future terms will also be zero; hence the summation 
may cease when the first zero term is obtained after a. + 1 steps. 

A final algorithm also requires a. + 1 steps and is essentially a variation 
of that if corollary 7 but requires the division of smaller integers and always 
divides by p^ rather than p ^ . 

8. Corollary. Let 


0(a) 



P i 

0i(a) 



P i 

02(a) 


p i 

03(a) 


represent a continued division algorithm. Then a. = k - 1 such that 0^(a) is 
the first quotient that is not an integer; at this point the algorithm terminates. 
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Clearly all of the above algorithms are predicated on our ability to 

k 

determine whether or not 0(a)/p is an integer. As such 0(a) must be stored 
with sufficient accuracy to determine this fact if a is to be recovered exactly. 
Now any integral error in 0(a) certainly results in the change of some prime 
factors; hence the integer 0(a) must be known exactly. Of course since it is 
known a priori that 0(a) is an integer, one can accept fractional errors of less 
than half, but no integral errors are acceptable if a is to be recovered exactly. 
On a computer, therefore, 0(a) must be stored as an integer number; it does 
not suffice to store a few significant figures together with an exponent. The 
application of the data compression scheme is thus limited by the size of the 
integers 0(a) that can be stored on the available computer. Of course multiple 
precision words can be used but even then the number 0(a) may become too 
large for the machine unless the vector a is in some way restricted. In the 
particular case of the vectors resulting from topological matrices a is sparse; 
hence many of the factors in 0(a) are unity and are composed of small inte- 
gers; hence the prime factors of 0(a) are raised to small powers only. These 
two conditions thus tend to keep 0(a) within reasonable bounds and render the 
procedure practicable for such matrices. Note that although the sparseness of 
the topological matrices is necessary to render the data compression scheme 
feasible, one does not need to assume a specified pattern of sparseness and 
may therefore manipulate compressed topological matrices without regard for 
variations in the sparseness pattern so long as the magnitude of 0(a) remains 
within reasonable bounds. 


Conclusions 

To implement the preceding data compression scheme one needs for the 
encoding process an ordered table of primes and for the decoding process a 

subroutine for evaluating j^l - 0(a) /p^ + [^0(a)/p^] J (that is an algorithm for 

k 

determining whether or not 0(a)/p is an integer) . In the former case one 
can, of course, simply tabulate a sufficiently long list of primes, but a more 
efficient approach is to use a polynomial approximation. That is a polynomial 
P(x) such that 

p. < P(i) < p. + 1 (62) 

1 i 

which yields p. exactly upon taking the integer part of P(i) . Since the primes 
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are a reasonably smooth monotonic function of the integers, a relatively low 
order polynomial P(x) can be used to calculate, exactly, a large number of 

primes. In the case of the decoder one can, of course, calculate 0(a) /p 
exactly prior to determining whether or not it is an integer, but since one is not 

ic 

really interested in the value of 0(a)/p , a simpler subroutine is possible 

ki 

wherein the calculation of 0(a) /p is carried only far enough to determine 
whether or not it is an integer. 
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